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O ; ABSTRACT 

1^ ■ Using the deep Spitzer Infrared Array Camera (IRAC) observations of the 

i/-^ ■ Great Observatories Origins Deep Survey (GOODS), we study the stellar masses 

^ ' and star formation histories of galaxies at ^ ~ 6. The IRAC instrument provides 

> ■ the best opportunity to estimate the stellar masses of galaxies at these redshifts 

J^ ■ because it samples their rest-frame optical fluxes, which are less prone to dust 

•/^ ' extinction and are more sensitive to the light from longer-lived stars. Our study 

(^ ■ is based on the zyrs-band dropout sample selected from the GOODS southern and 

^ ■ northern fields (~ 330 arcmin^ in total), several of which already have spectro- 

■"^ ' scopic confirmations. In total, we derive stellar masses for 53 i775-band dropouts 

P^' that have robust IRAC detections. These galaxies have typical stellar masses 

O ■ of ~ 1O^°M0 and typical ages of a couple of hundred million years, consistent 

-t— » ' with earlier results based on a smaller sample of 2; ~ 6 galaxies in the Hubble 

ci • Ultra Deep Field. We also study 79 zyys-band dropouts that are invisible in the 

<^ • IRAC data and find that they are typically less massive by a factor of ten. These 

^ ■ galaxies are much bluer than those detected by the IRAC, indicating that their 

j^ ■ luminosities are dominated by stellar populations with ages < 40 million years. 

We discuss various sources of uncertainty in the mass estimates, and find that 
our results are rather robust. The existence of galaxies as massive as 10^^ Mq at 
z ^ Q can be explained by at least one set of N-body simulations of the hierarchi- 
cal paradigm. Based on our mass estimates, we derive a lower limit to the global 
stellar mass density at z ^ 6. Considering the range of systematic uncertainties 
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in the derived stellar masses, this lower limit is 1.1 to 6.7x10^Mq Mpc^^ (co- 
moving), which is 0.2 to 1.1% of the present-day value. The prospect of detecting 
the progenitors of the most massive galaxies at yet higher redshifts is explored: 
a deep, wide-field near-IR survey using our current technology could possibly 
result in positive detections at z > 7. We also investigate the implication of our 
results for reionzation, and find that the progenitors of the galaxies comparable 
to those in our sample, even in the most optimized (probably unrealistic) sce- 
nario, cannot sustain the reionization for a period longer than ~ 2 million years. 
Thus most of the photons required for reionization must have been provided by 
other sources, such as the progenitors of the dwarf galaxies that are far below 
our current detection capability. 

Subject headings: cosmology: observations — galaxies: evolution — galaxies: 
luminosity function, mass function — infrared: galaxies 



1. Introduction 

The Infrared Array Camera (IRAC; Fazio et al. 2004) of the Spitzer Space Telescope 
(Werner et al. 2004) has opened up a new window for the study of galaxies at very high 
redshifts. Its impressive sensitivity in the 3.6 and 4.5 /xm channels enables the detection 
of galaxies as distant as z ~ 6, and thus, for the first time, offers the unique opportunity 
of investigating the rest-frame optical properties of galaxies in the early universe. At these 
wavelengths, the light from a galaxy is less affected by dust extinction and is more sensitive to 
the longer-lived stars that dominate the stellar mass. For these reasons, IRAC observations 
can provide a wealth of information about the stellar population of galaxies at very high 
redshifts, such as their ages and stellar masses. 

The past year has witnessed substantial progress in this frontier. Egami et al. (2005) 
used IRAC to detect a galaxy at z ^ 6.7 lensed by a foreground cluster and estimated that 
this object has a stellar mass of ~ IO^Mq and is at least ~ 50 Myr old. Eyles et al. (2005) 
discussed two field galaxies at z ^ 6 that are detected in the Great Observatories Origin 
Deep Survey (GOODS) IRAC data and found that they are a few hundred Myr old and 
their stellar masses are a few xlO^°M0. Yan et al. (2005; Paper I) analyzed three galaxies at 
^ Ri 6 and 11 galaxies at z ^ 5 in the Hubble Ultra Deep Field (HUDF) that are detected in 
the GOODS IRAC observations, and concluded that the observed number density of massive 
galaxies (> 10^*^ M0) is consistent with predictions from contemporary numerical simulations 
of galaxy formation in a cold dark matter (ACDM) universe. Chary, Stern & Eisenhardt 
(2005) investigated the red IRAC color of a Lya emitter at z = 6.56, and suggested that 
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the red color could be caused by the presence of a very strong Ha emission line (rest-frame 
equivalent width ~ 0.2yum) redshifted to the 4.5/im channel. Mobasher et al. (2005) studied 
a galaxy in the HUDF that is not visible at optical wavelengths but is very prominent in the 
IRAC images, and suggested that it could be an old galaxy at z ^ 6.5 with a surprisingly 
high stellar mass of 5.7 x lO^^M© (but see Chen & Marzke (2004) and Yan et al. (2004) for 
different interpretations of this same object). 

While some of the conclusions from these studies are still tentative, they show that deep 
IRAC observations have enabled a new perspective on studies of early galaxy formation. 
Previous studies, however, have been limited by the small number of galaxies that were 
analyzed, and by the small fields of view, where clustering variance might substantially 
affect space densities or other properties of the galaxy samples. In this paper, we extend the 
analysis to the combined areas of the two GOODS fields, in the regions of the Hubble Deep 
Field North (HDF-N) and Chandra Deep Field South (CDF-S). These fields cover roughly 
30 times more solid angle than the HUDF alone, along two sightlines. This substantially 
increases the number of galaxies we can study, and should help average over variance due to 
large scale structure. We select F775W-band dropouts (hereafter, iyys-dropouts) from the 
i/S'T Advanced Camera for Surveys (ACS) GOODS observations (Giavalisco et al. 2004a), 
and identify counterparts from the full-depth, two-epoch Spitzer IRAC observations of these 
fields (Dickinson et al., in preparation). These candidate z ~ 6 galaxies are presented in 
§2 and their IRAC properties are discussed in §3. We estimate their stellar masses in §4 
and discuss the implication of our results in §5. Our conclusions are summarized in §6. 
The magnitudes quoted in this paper are all in the AB system. Throughout the paper, we 
adopt the following cosmological parameters based on the first-year Wilkinson Microwave 
Anisotropy Probe (WMAP) resuhs from Spergel et al. (2003): ^m = 0.27, I^a = 0.73, and 
Hq = 71 km s~^ Mpc~^ (our results are not affected if using the values from the three-year 
WMAP results for these parameters). All volumes quoted are co-moving. 



2. Galaxy Candidates at z ?5i 6 in the GOODS Fields 

The z ^ 6 galaxy candidates used in this paper are the irys-dropouts selected from the 
complete, 5-epoch GOODS ACS observations in both the HDF-N and CDF-S (Giavalisco et 
al. 2004a). The selection of these candidates will be discussed in detail by Giavalisco et al. 
(2006, in preparation). Briefly, the selection is aimed at the redshift range of 5.5 ^ z ^ 6.5, 
and requires a candidate have a red color across the F750W (^775) and F850LP (zs^o) bands, 
and be invisible in the F606W (V^eoe) and F435W (-B435) bands. 

Quantitatively, the candidates must meet these criteria: S/N > 5 in zgso, (^775 — -^sso) > 
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1.3 mag, and S/N < 2 in both the Veoe and -B435. These criteria are the same as have 
been previously used by the GOODS team (Giavahsco et al. 2004b; Dickinson et al. 2004). 
While the colors are calculated based on the SExtractor (Bertin & Arount 1996) MAGJSO 
magnitudes, MAG_AUTO magnitudes are used when quoting the Zgso magnitudes of these 
candidates as they are closer to the "total" magnitudes. As brown dwarfs in our Galaxy have 
similar red colors, an "image sharpness" criterion, based on comparing peak zgso-band surface 
brightness within a 0'.'125 aperture to the isophotal magnitude, was used to reject possible 
point sources from the iyys-dropout sample down to 2:850 = 26.5 mag. This criterion was 
purposefully tuned to be conservative: it may improperly reject a few galaxies, but should 
successfully eliminate most stars down to that magnitude limit. At fainter magnitudes, it 
becomes difficult to robustly distinguish stars and compact galaxies in the GOODS ACS 
images, as the stellarity loci of stars and galaxies start to merge for any of the criteria we 
have considered. 

After visual inspection to eliminate likely spurious sources, the resulting sample consists 
of 274 i775-dropout candidates (145 in the south, 129 in the north). Of these, eight more 
were further rejected as possible stars on the basis of having full-width-at-half-maximum 
(FWHM) < 0'.'125; two of these have been spectroscopically confirmed as stars by Vanzella 
et al. (2006). After this process, 142 and 124 valid candidates remain in the CDF-S and 
HDF-N, respectively. 

We have identified five close pairs among these iy/s-dropouts, three in the CDF-S sample 
and two in the HDF-N sample, respectively. The separation of the components in these pairs 
are all smaller than 2", which means that these pairs, if detected in the IRAC images, will 
be detected only as single sources (see below). Therefore, for the statistical purpose of 
this paper, we take the number of iyys-dropouts in the CDF-S and HDF-N as 139 and 122, 
respectively. 



3. IRAC Detections of z ~ 6 Galaxy Candidates 

The IRAC data used in this paper are the mosaics of the full, two-epoch data of the 
GOODS observations in both the CDF-S and HDF-N (Dickinson et al., in preparation; see 
also http://data.spitzer.caltech.edu/popular/goods/Documents/goods_dataproducts.html). The 
mosaics cover approximately 10 x 16 in each field (total coverage of ~ 330 arcmin^), and 
have a nominal exposure time of ~ 23.2 hours in each of the 3.6, 4.5, 5.8, and 8.0 /xm chan- 
nels. Because of the way the observations were designed, the exposure time in the central ~ 
20% stripe of each field is twice as long. The final pixel size of the mosaics is 0.6 pixel"^, 
i.e., roughly half of the native IRAC pixel size. Sources were detected in a weighted sum 
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of the 3.6 and 4.5 /iin images using SExtractor, and their magnitudes were measured in 
each band through circular apertures that are 3.0 in diameter. These aperture magnitudes 
were converted to total magnitudes after adding aperture corrections (see Paper I). We only 
include sources with S/N > 3 (measured in 3.6/im within the above mentioned aperture) for 
cross-matching. 

The i775-dropouts were cross-matched with the IRAC sources following the procedures 
adopted in Yan et al. (2004; 2005). As the FWHM of the point spread function (PSF) in 
the 3.6/im channel is about 1.8", the matching was done with a 2" search radius. This radius 
was chosen partly because of the close pairs in our sample, whose IRAC centroids are offset 
from those of the individual components. The matched sources were then visually inspected 
to ensure that the identifications were secure. In order to avoid ambiguity in interpreting the 
measured fluxes, we excluded any IRAC sources that are blended with neighbouring objects. 

In total, 64 irys-dropouts have thus been securely identified in the IRAC data, with 40 
in the CDF-S and the remaining 24 in the HDF-N. Among these sources, there are two pairs 
in the CDF-S and one pair in the HDF-N, respectively, all of which are counted as single 
IRAC objects in the numbers quoted above. Except for these three pairs, all other objects 
have their IRAC and ACS centroids matched to within 0.6 , which is fully consistent with 
our previous experience (see Yan et al. 2004; 2005). 

The Lyman-break selection of galaxies at 2; ^ 6 is known to suffer some contamina- 
tion from low-redshift, early- type galaxies that have similar red colors. Our IRAC data 
can be used to further reduce these contaminators. A number of the iyys-dropouts iden- 
tified above have very bright counterparts in the IRAC images, giving flux density ratios 
fu(y3-GfJ''>TT')/fu{zs5o) > 20 (or equivalently, (2:850 — ms^^rn) > 3.25 mag). Thus they satisfy 
the selection criterion of the "IRAC-selected Extremely Red Objects" (lEROs; Yan et al. 
2004), which are likely old, passively-evolving galaxies at z ^ 2.4. Therefore, these objects 
are excluded from our final sample. There are six and five such lEROs among the IRAC- 
detected iyrs-dropouts in the CDF-S and HDF-N, respectively. The ratio of contamination 
due to lEROs is 15% for the CDF-S (six out of 40) and 21% for the HDF-N (five out of 24). 

We also consider the iyys-dropouts that are not detected by the GOODS IRAC obser- 
vations. We visually examined the IRAC images at the positions of the i775-dropouts that 
were not cross-matched in the above procedure and found that 45 objects in the CDF-S and 
34 objects in the HDF-N were undetected because of their faintness — quantitatively, they 
all have S/N < 3 in both the 3.6 and 4.5 /im images. The remaining 118 objects were not 
cross-matched because they are blended with nearby, unrelated sources. 

To summarize, we are left with a final sample of 250 objects for analysis (133 in the 
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CDF-S and 117 in the HDF-N). We have counted the close pairs as single sources, and have 
excluded the possible contaminators (brown dwarfs and lEROs). The final sample of 2775- 
dropouts that are securely detected in the IRAC data has a total of 53 objects (34 in the 
CDF-S and 19 in the HDF-N). The 3.6 /im-band magnitudes of these objects range from 
23.27 to 26.50 mag. The final sample of i775-dropouts undetected in the IRAC data consists 
of 79 objects (45 in the CDF-S and 34 in the HDF-N). In the rest of the paper, we will 
refer to these two samples as the "IRAC-detected sample" and the "IRAC-invisible sample" . 
The remaining 118 i775-dropouts (54 in the CDF-S and 64 in the HDF-N) are blended with 
unrelated neighbours. lEROs make up ~ 17% of the non-blended sample, and it is plausible 
that the same fraction also holds in the blended sample. Therefore, a total of 98 objects 
in the blended sample could be galaxies at z ^ 6 after excluding this fraction of lERO 
contamination. We will refer to the sample of these 98 objects as the "blended sample" . 

Despite that we have utilized every means to eliminate various sources of contamination, 
our final sample could still have some contaminators. This is particularly possible among 
the objects of low S/N, where spurious sources are more likely to occur and low-redshift, red 
galaxies (not necessarily as extreme as the lEROs) are more easily scattered above the color 
selection threshold. To address this issue, we use the much deeper ACS data in the HUDF 
to assess the quality of our candidate selection. Thirteen objects in our CDF-S candidate 
list are within the coverage of the HUDF, and all of them are verified to be real objects. 
Therefore, the contamination caused by spurious sources seems to be negligible to the S/N 
level that we consider. However, three of these thirteen sources are not among the HUDF 
i775-dropout list of Yan & Windhorst (2004b). These three sources {S/N ^ 5 in the GOODS 
data) have ^775 — -Zsso < 0.5 mag in the deeper HUDF data, indicating that the contamination 
rate in our final sample caused by photometric error could be as high as ~ 23%. We will 
take this into account when discussing the implications of our results in §5. 

Fig. 1 shows the Zgso-band magnitude distributions of the IRAC-detected and IRAC- 
invisible samples as histograms of different colors. Their median magnitudes are 26.32 and 
27.00 mag, respectively. For comparison, the median magnitude of the blended sample is 
26.71 mag. While the objects in the IRAC-invisible sample are fainter on average, the 
majority of them are still within the range of the IRAC-detected sample. 



4. Constraining the Stellar Masses of z ~ 6 Galaxies 

In this section, we constrain the stellar masses of z ~ 6 galaxies in both the IRAC- 
detected and the IRAC-invisible samples as defined in §3. For this purpose, we compare our 
observations to the stellar population synthesis models of Bruzual & Chariot (2003; hereafter 
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BC03). As per Paper I, we explore the models of the following star formation histories 
(SFHs): instantaneous bursts (or Simple Stellar Populations, SSPs), and continuous bursts 
with exponentially declining star formation rate (SFR) in the form of SFR ex r~^ exp{t/T), 
where the time-scale r ranges from r = 10 Myr to 1 Gyr. The step size in r is 10 Myr 
for r = 10 Myr to 0.1 Gyr, and is 0.1 Gyr for r = 0.1 to 1 Gyr. For each of these SFHs, 
we generate models with ages ranging from 1 Myr to 1 Gyr; the step size in age is 10 Myr 
when the age is less than 0.1 Gyr, and is 0.1 Gyr when it ranges from 0.1 to 1 Gyr. To 
facilitate comparison to other, published studies, we use a Salpeter initial mass function 
(IMF; Salpeter 1955) extending from 0.1 to 100 Mq. Many recent studies indicate that 
the actual IMF in most environments probably has fewer low-mass stars (< IMq) than are 
predicted by the Salpeter IMF power-law slope (e.g., Kroupa 2001; Chabrier 2003). These 
stars contribute very little to the observed light from galaxies, and therefore, to first order, 
changing the low-mass IMF simply rescales all masses and star formation rates derived from 
comparison to stellar population models. Generally speaking, using a Chabrier IMF would 
lower the mass estimates by about 50-60%. Changes to the IMF shape and slope at higher 
masses, however, may lead to wavelength-dependent and time-dependent changes in the 
mass-to-light ratios (M/L) relative to the assumed Salpeter models. 

In Paper I, we carried out detailed SED analysis using models of different metallicities 
and various reddening values, and allowed the redshift be a free parameter when it was un- 
available. When single-component models did not fit well, we also considered two-component 
models. Here, however, it is difficult to take this approach, because the objects in the IRAC- 
detected sample are only significantly detected in two to three passbands (-2350 and 3.6/im, 
some are also detected in 4.5/im). This is particularly true for the IRAC-invisible objects, as 
they are only significantly detected in one band (-2350) • Therefore, we have to use a different 
method, detailed below. We will mainly concentrate on models of solar metallicity and zero 
reddening, but will also discuss the effects of different metallicities and reddening values. 

Currently, seven objects in the IRAC-detected and -invisible samples have already been 
spectroscopically confirmed aX z > 5.5 (Stanway et al. 2004a, 2004b; Dickinson et al. 2004; 
Vanzella et al. 2006; Stern et al., in preparation; Vanzella et al., in preparation). For the 
sake of simplicity, however, we choose to assign a common redshift oi z = 6 for all objects. 
We will discuss the bias caused by this simplification in §4.3 together with other sources of 
uncertainty. 



4.1. Stellar Masses of IRAC-Detected Objects 

For each object, we first examine how its luminosity in the IRAC bands constrains the 
upper limit of its stellar mass. The luminosity of an object with fixed mass is determined 
by its age and SFH, and becomes fainter as it gets older; i.e., its M/ L becomes larger as it 
evolves. It can reach the largest possible M/ L if it forms all its stars through an instantaneous 
burst and passively evolves to its maximal age (i.e., a SSP; see Fig. 2). Therefore, an 
object can have the highest possible mass (assuming no dust extinction) if its luminosity is 
dominated by such a "maximally old", single-burst component. This methodology is similar 
to that adopted in Papovich et al. (2001). In practice, we assume that all the fiux of 
an object detected in the 3.6/im channel comes from this maximally old component. The 
3.6/im channel is chosen because it has the highest sensitivity of the IRAC bands, and thus 
the measurements in this channel have the highest S/N. As we assume the objects in our 
sample are all at z = 6, their maximal age allowed (i.e., setting the formation redshift 
Zf = oo) in our adopted cosmology is 0.95 Gyr. To match the model grid in age domain, we 
take their maximal ages as 1.0 Gyr. The stellar mass upper limit of a given object is then 
obtained by comparing its 3.6/im magnitude to the prediction given by the 1.0 Gyr-old SSP 
model. Hereafter we refer to this limit as "Mmax"j and the method as the "maximally-old 
component method" ^. 

Next, we find the lower bound to the stellar mass (hereafter M^[^) of each object. If the 
SFH of an object is known, its (-2850 —""^3.6^™) color can be used as an indicator of its age (see 
Fig. 3). Its mass can be derived through comparison of the observed 3.6/im magnitude to 
the prediction of the corresponding model at the inferred age. As we do not know what SFH 
a given object has, we consider the full range of models and obtain one mass estimate for 
each model. For each galaxy, the minimum among all these estimates is taken as its Mmin- 

Given our assumptions, the true mass of a galaxy will lie between Mmin and Mmax- 
While there is no other information that can further help us judge what its true mass might 
be, we use the median mass from the set of exponential SFH models fitted to match the 



^Assuming no dust extinction (discussed in §4.3), the Mmax value thus obtained marks a firm upper limit 
to the mass. The maximally-old component always underpredicts the flux in zsso-band, and to explain the 
observed zgso-band flux would require a secondary, young component be added. In such a two-component, 
"old-|-young" model, the derived mass of a galaxy will be smaller than the Mmax value that we obtained, 
because the contribution from the young component will reduce the M/L. The exact amount of such a 
reduction depends on the assumed age and SFH of the young component; using the suite of models shown 
in Fig. 3 for this young component, the derived maximum mass will be smaller than our current Mmax value 
by a few to 50% on average. As wc have no further constraints on the young component, we choose to adopt 
the Mmax value based on our current simplistic model as a safe upper limit. 
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(-2850 ~ iTT'S.Gfim) color. We refer to this as the "representative mass" of this galaxy (hereafter 
^rep)- We also refer to the age corresponding to this Mrep as its "representative age" (T^ep)- 
The representative masses range from 0.09-7.0 xlO^^M© (the median is 9.5x10^Mq) and the 
representative ages range from 50-400 Myr (the median is 290 Myr); both are consistent with 
one of the general conclusions in Paper I that some galaxies as massive as a few xlO^^M© 
were already in place at z ~ 6, and that they are typically a few hundred million years old. 
The histograms of the three values, Mmin, Mrep and Mmax, are shown in the top panel of Fig. 
4. As a further justification of using M^^p and T^p, we note that Mj-ep for J033240. 01-274815.0 
is 2.1xlO^'^M0, which is well within a factor of two of the best-fit mass of 3.4xlO"'^°M0 that 
Paper I derived for this object using a more sophisticated SED analysis. For comparison, 
the Trep value of this object is 200 Myr, while the analysis of Paper I shows that the evolved 
component of this object has an age of 500 Myr. 



4.2. Stellar Masses of IRAC-Invisible Objects 

The fact that many objects are undetected in the IRAC data immediately suggests that 
they are likely much less massive that the IRAC-detected ones. The 3.6jum flux density upper 
limits of these sources, measured as 2 o" fluctuation within a 3 -diameter aperture, range from 
0.029 to 0.091 /iJy (27.73 to 26.50 mag), with the median of 0.042 /xJy (27.35 mag). Using the 
maximally-old-component method described above, we derive the stellar mass upper limits 
of these objects based on their flux density upper limits. Again, we assume all galaxies are at 
z = 6. These upper limits, which are shown in the dotted histogram displayed in the lower 
panel of Fig. 4, range from 3.0-9. 4x10^Mq. The median of these upper limits is 4.3x10^ Mq, 
which is not only less than the upper limit of the IRAC-detected sample in general, but is 
also signiflcantly less than the median of the representative mass of the latter by a factor of 
two. 

Fig. 5 compares the {zg^Q — m^Q^rn) colors of the IRAC-detected objects (red squares) 
with the upper limits of these colors of the IRAC-invisible objects (downward arrows)®. 
This comparison suggests that the IRAC-invisible objects are not only less massive than the 
IRAC-detected ones but are also bluer on average, indicating that young stellar populations 
are playing a dominant role — by luminosity — in these galaxies. 



®The "gap" in this figure that seemingly seperates the IRAC-detected and IRAC-invisible objects does 
not imply that these are two distinctly different populations. Instead, the transition between these two 
samples should be rather smooth. The "gap" is largely caused by the sources that fall just slightly below 
our detection threshold in IRAC {S/N > 3), whose (zsso — J^ia.g/^m) color limits are then calculated by using 
their 2 a flux upper limits in the 3.6/im-band. 
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To further test the robustness of the conclusion that these objects are less massive, 
we stacked their 3.6/im images to increase the S/N of the detection. For each object, a 
12.6 X 12.6 section (i.e., 21-pixel on a side) around its center was cut out from the 3.6/im 
image. We then combined these image "stamps" together by taking their median value at 
each pixel. We adopted median rather than average as the stacking algorithm, because the 
former is more effective in rejecting the contaminating pixels from the neighbours around 
invisible sources. The final median stack is displayed in the left panel of Fig. 6, which shows 
a visible, albeit weak, source at the center (see pixels within the circle). For comparison, 
the right panel of this figure shows the stack of the same number of randomly chosen image 
stamps, where no detectable source at the center can be seen. 

This weak source represents the average property of the IRAC-invisible iyys-dropouts 
in the 3.6/im channel. To extract its flux, we utilized a 1.2"-diameter (2 pixels) aperture. 
The background value and the aperture correction to total flux were determined through 
simulation. The simulation was performed for 20 runs, and in each run 830 different artificial 
galaxies with a total magnitude of 27.50 mag were randomly distributed over the real image. 
Image stamps of these artificial galaxies (16,600 in total) were cut out and stacked in the 
same way as we stacked the IRAC-invisible sources. The peak value of the pixel histogram 
of this stack was chosen as the background value, and the magnitude of the stacked artificial 
galaxy was extracted using a 1.2"-diameter aperture. The aperture correction was calculated 
as the difference between this aperture magnitude and the input total magnitude. The same 
background value and the aperture correction were then applied to the real stack of the 
IRAC-invisible object. Its total magnitude thus obtained is 27.44 mag, which is consistent 
with the flux density upper limits described above. 

By the same token, the representative property of this median IRAC source in Zs^o-hand 
can be described by the median 2:850 magnitude of the IRAC-invisible sample, which is 27.00 
mag. Using the same analysis in §4.1, we obtain (Mmin, M-ep, ^max) estimates for a typical 
source in the IRAC-invisible sample as (1.5 x 10^, 2.0 x 10^, 5.9 x 1O^)M0. For all three types 
of estimates, the masses derived from the median-stack photometry are much smaller than 
those for the objects in the IRAC-detected sample (see the dashed lines in the top panel of 
Fig. 4). The T^ep value of this median source is 30 Myr. 



4.3. Sources of Uncertainty in Mass Estimates 

A major source of uncertainty is in the systematic of the measurement of the (-2350 — 
fns.Gfmi) color. As it is used as the age estimator, this color affects both M^in and Mj-cp (but 
not Minax)- In this study, we use total magnitude for m^Q^rn and SExtractor MAG_AUTO 



-11- 



for ^850- Although MAG_AUTO is frequently taken to represent the total magnitude of a 
galaxy, in practice it may be subject to biases, particularly for faint objects in high-resolution 
HST images. Monte Carlo simulations of photometry for artificial objects in the GOODS 
ACS images indicates that at the typical magnitudes of our iyys-dropouts, the MAG_AUTO 
magnitude underestimates the total ^rgso-band flux of an object by ~ 0.5 mag. This means 
that the true (zsso — "n^s.e^m) color could be 0.5 mag bluer than our current estimate, and 
thus the age of a given object could be younger, and hence the M/ L could be smaller. We 
find that our current M^^^ and M^ep values could have been overestimated by ~ 21% if the 
the (2:850 — rriss^m) color is indeed 0.5 mag bluer than we adopted. 

In addition to photometry, there are more complicated sources of uncertainty. When 
deriving the results shown in previous sections, we have made three major simplifications, 
namely, the galaxies (or equivalently, the models) are all at z = 6, all have solar abundance 
and all are free of dust obscuration. Here we discuss the systematics introduced by such 
simplifications. We mainly discuss the effects on the derived stellar mass, as our current 
study largely concentrates on this quantity. We consider how the derived masses would 
change if we vary the assumed values for the redshift, metallicity, and reddening. We express 
the systematics in terms of relative difference in mass, AM^/M, where the superscript i can 
be "2;" (redshift), "Fe/H" (metallicity) and "red" (reddening). The quantity AM is defined 
as AM = M — M , where M is the value obtained using our adopted models, and M is 
the value obtained when the model parameter in question is changed. A positive value of 
AM^/M indicates that our default models yield a larger estimate for the mass than those 
where parameter i is changed, and a negative value means the opposite. We consider the 
effects of model parameter changes on all three mass estimates considered here, namely, 
Mmin, Mrep, and Mmax- We ouly discuss the IRAC-detected sample, but the conclusions 
could similarly be applied to the IRAC-invisible sample. 

By fixing the redshift, we ignore the differences in the redshifted SED, in the luminosity 
distance, and in the amount of the IGM absorption. The first two factors (especially the 
luminosity distance) are the most relevant in deriving Mmax using a maximally-old SSP. 
Using a z = 6 SSP will overestimate the Mmax values for objects at lower redshifts, and 
will underestimate them for objects at higher redshifts. This is demonstrated in the bottom 
panel of Fig. 7 for two extreme cases oi z = 5.5 (asterisks) and z = 6.5 (open squares). 
The Mmax values of our objects in the IRAC-detected sample have been recalculated using 
a SSP at z = 5.5 and a SSP at z = 6.5, respectively. Our original values would overestimate 
by 21% if the objects were actually all at z = 5.5, and would underestimate by 27% if they 
were actually all at z = 6.5. 

For galaxies at z > 5.9, the Lya forest suppresses the fiux detected in the Z85o-band. 
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This has no effect on the maximal mass estimates, Mmax, which is based only on the ffux in 
the 3.6/xm channel. However, it does impact the mass estimates Mmin and Mj-ep, which make 
use of M/ L derived from the Zgso-S.G/xm color. For example, if the object under question is 
actually at 2; > 6, its observed (2:850 — ""^3.6^™) color would be redder than what is expected 
when at z = 6, because its flux in the Zgso-band is more severely suppressed. Applying a 
z = Q template to such a color will result in an artificially older age (see Fig. 3), and hence an 
artificially larger M/ L (see Fig. 2). As a result, both M^i^ and M^p will be overestimated. 
This overestimate, however, is partially canceled by the change in the luminosity distance 
and the /^-correction. The top and middle panels of Fig. 7 show that the assumption of 
z = Q would overestimate Mmin and Mj-cp by about 25% for galaxies that were in fact at 
z = 6.5 (see the open squares). At 2; < 6, this effect results in underestimates. However, as 
the change of the IGM opacity across the zgso-band is negligible at this redshift regime, the 
mass estimates change by only a few percent if the galaxies were instead aX z = 5.5. This is 
shown in the top and middle panels of Fig. 7 as asterisks. 

In Paper I, we considered the metallicities varying from solar abundance to 1/200 of 
solar, the latter being the smallest value available from the BC03 models. We found that 
the z ^ Q objects in the HUDF are usually best fitted by models with solar abundance. For 
this reason, we have adopted the simplifying assumption of solar metallicity for our fiducial 
models. As it is not likely that objects in the early universe would have higher abundance, 
here we only examine the extreme case at the low metallicity end, i.e., how our results would 
be biased if all the objects actually have their abundance of 1/200 of solar. Fig. 8 shows 
the relative differences of the three mass estimates, all evaluated at 2; = 6. Models with 
subsolar metallicites generally have bluer colors at a fixed age; hence, older ages are required 
to match the observed colors, and this leads to larger M/L. Therefore, the masses based on 
the colors (Mmm and Mj-ep) derived using the solar abundance models are smaller than those 
derived from models with Z = Zq/ 200 by 70% and 32%, on average. For Mmax, instead, the 
solar abundance models have higher M/L, and hence the masses are 32% larger than those 
for models with 1 /200th solar metallicity. 

Dust reddening and obscuration will also affect the derived masses. In Paper I, we 
found that the best-fit models always have near-zero reddening, although an amount up to 
E{B — V) = 0.4 to 0.5 mag could not be excluded. This is because we had to account for the 
very blue rest-frame UV colors observed for the galaxies studied in that paper, which we did 
by invoking two-component models with both younger and older stars. The best-fit models 
always have E{B — V) ^ if we require both components to have the same reddening. 
However, if this ad hoc requirement is removed, the reddening becomes unconstrained. Here, 
we consider the effects of non-zero reddening on our results. For mass estimates based on 
a color measurement, i.e., Mmin and Mrep, the effects of reddening and age on M/L derived 
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from the models tend to cancel, resulting in relatively small changes in the final derived 
masses. Dust makes the observed 3.6/im flux fainter, thus increasing M/ L. However, it 
also reddens the spectrum, requiring a younger model to match the observed (zgso — 3.6/im) 
color, and this younger model will have an intrinsically smaller M /L (before applying the 
extinction). For illustration, we consider the effects of E{B — V) = 0.2 mag of extinction, 
assuming the attenuation law of Calzetti et al. (2000). The results are shown in Fig. 9. 
Although this reddening corresponds to nearly 1 mag of extinction in the observed 3.6//m 
band and more than 2 mag at Zgso, the average change in the derived Mmin and Mrep is only 
13% and 30% (overestimated), respectively, as the model age changes to compensate the 
extinction. Mmax, on the other hand, is derived from the photometry at 3.6/im alone, and 
hence changes by a factor of 3.3 (underestimated). 

In the end, we flnd that most of our mass estimates are relatively robust against large 
changes in the values assumed for the redshift, metallicity, and extinction. In most cases, 
varying those parameters over the ranges described here, the derived masses change by less 
than 35%. If metallicity were as low as l/200th solar, then our values for Mjam (derived 
assuming Zq models) would be underestimated by 70% on average. The largest effect would 
be that of reddening on our estimates of Mmax? which could be increased by a large factor if 
reddening were substantial. However, we note that this would imply substantial reddening 
on the old (1 Gyr) stellar population assumed in the toy model from which the mass upper 
limits are derived. In practice, the mass could be increased arbitrarily by adding dust to 
an old stellar population, but such a model seems rather arbitrary and unphysical, and also 
fails to take into account the observed UV light and UV-to-optical rest-frame colors of the 
galaxies. Instead, the mass estimates derived from the observed colors (Mmin and Mrep) are 
comparatively insensitive to reddening. Other stellar population properties are more sensitive 
to these systematic effects. In particular, both age and SFR are completely degenerate with 
extinction when only a single color is used for analysis. Thus, while we have shown that 
the mass estimates are comparatively robust, other properties are more uncertain, and rely 
more heavily on the inference for low extinction from more detailed analysis of the smaller 
HUDF sample in Paper I. 



4.4. Legitimacy of Adopted Declining SFHs 

The stellar mass estimates that we obtained in this study hinge upon the exponentially 
declining SFHs that we assumed. Such SFHs imply that the star-formation activities of these 
objects were more intense in the past than in the epoch when they are observed. While it 
is impossible to directly verify this picture based on our current data, we can at least test 
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if our results are self-consistent with the assumption that the SFHs of these objects are a 
declining function of time. 

One quantity that can be used to address this question is "specific SFR" , which is defined 
as the ratio of the current SFR and the total stellar mass. The reciprocal of this quantity is 
the period of time (denoted as T) during which a galaxy could build up its total mass if it 
were maintaining the same, constant SFR in the past. The current SFR can be calculated 
according to eqn. 2 of Madau, Pozzetti & Dickinson (1998), i.e., Luv = 8-0 x 10^'^ x SFR, 
where Luv is the flux (in ergs/s/Hz) at rest-frame 1500A. The SFR thus derived is insensitive 
to the past SFH when the age of the galaxy is much longer than the life time of O and B 
stars at the main sequence. We convert the Zgso-band magnitudes to Luy, again assuming 
that the objects are all at z = 6, and ignoring the caveat that the Zgso-band actually samples 
the continuum at around 1300A. 

Fig. 10 shows the specific SFR (calculated by dividing the Mj-cp into the SFR value) 
of the IRAC-detected sample as the function of stellar mass (solid triangles). The result 
for the median-stack of the IRAC-invisible objects is also plotted (the asterisk). The data 
points form a tight sequence that is between constant SFR of 3.0 and 30 MQyr^^, which are 
shown as the two straight lines to the left and right, respectively. This primarily reflects the 
rest-frame UV magnitude range of the sample (24.9 < zg^o < 27.4 mag) and its conversion 
into star formation rates (30 > SFR/{MQjr~^) > 3). The lower envelope in specific star 
formation rates corresponds to the ^§50 niagnitude limit of the sample. The upper envelope 
may have greater physical significance, corresponding to the high-luminosity cut-off in the 
UV luminosity function for galaxies at z ~ 6. 

The T values are much larger than the representative galaxy ages, Trep, obtained in the 
previous section. More than half of these objects (58%) have T larger than the age of the 
universe at z ^ 6 (some are even > 4 Gyr), which is a clear indication that there would not 
have been sufficient time for them to acquire their masses if their past SFRs were as low 
as their current values. In other words, their SFRs must have been declining and must be 
much larger in the past, which is consistent with the assumed models. 

One caveat of the discussion above, however, is the possible role of dust extinction. If 
the galaxies in our sample suffered a reddening oi E{B — V) = 0.2 mag (as discussed in §4.3), 
their true ongoing SFR would be higher than the current estimates by a factor of ~ 8.5, and 
yet their M^^p values would be lower by ~ 30% (see §4.4). As a result, their T values would 
be smaller by a factor of ~ 12, and hence the majority of them could still build up their 
stellar masses within the time allowed by the age of the universe. However, even in this case, 
the T values (with the median of 155 Myr) are still much larger than Trep derived based on 
the reddened models (with the median of 10 Myr). Therefore, this is still consistent with 
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our assumption that the objects in our sample have higher SFRs in the past. 



5. Cosmological Implications 

The results presented above show that some galaxies as massive as a few xlO^^ Mq 
were already in place hj z ^ 6, and that their typical ages are a few hundred million years 
old, which are consistent with our earlier results in Paper I. In this section, we discuss the 
cosmological implications of these results. 



5.1. Comparison to ACDM Simulations 

Following Paper I, here we further investigate how well the observed number density of 
very massive galaxies aX z ^ 6 can be explained by the numerical simulations of the hierar- 
chical paradigm. As in Paper I, we make comparison to the mass functions (MF) predicted 
by two types of hydrodynamic simulations in a ACDM universe, namely, a Smoothed Parti- 
cle Hydrodynamics (SPH) simulation and a Total Variation Diminishing (TVD) simulations 
(Nagamine et al. 2004; Night et al. 2006). The simulation box sizes are 100h~^ Mpc and 
22h~^ Mpc, respectively. The cumulative MF from these two simulations are given in Fig. 
11. As the simulation volume is not unlimited, both MF are cut-off at the high- mass end (at 
around 1.8 x 10^^ and 2.0 x 10^° Mq for the SPH and TVD models, respectively). Neverthe- 
less, comparing to these models is still meaningful, as they extend to the regime where our 
observations are rather complete. For reasons as described in detail in Paper I, our IRAC 
data should be complete for detecting galaxies with M > 1.6 x IO^^Mq, which is indicated 
in the figure as the vertical dotted line. 

In Paper I, we derived a lower limit for the cumulative number density above the same 
mass threshold, which is 8.0xl0~^ Mpc~^, and we concluded that the predictions from 
the SPH and TVD simulations were a factor of 2.6 x and 4.6 x lager than this lower limit. 
However, we recently discovered that the mass functions from the SPH and TVD simulations 
were presented with incorrect normalization in Paper I, leading to predicted number densities 
that were too high by a factor of 6.0 and 3.8, respectively. Using the corrected numbers, the 
SPH and TVD predictions should be 0.4x and 1.2x of the observed lower limit. 

As the sample in Paper I is very small (three objects), the much larger sample in 
our current study offers a better statistics for comparison. Using the redshift selection 
function for our color criteria, which we evaluate from Monte Carlo simulations following 
the procedures described in Giavalisco et al. (2004b), we obtain the effective volume of 
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8.0x10^ Mpc'^ for our iyys-dropout sample. The cumulative number densities inferred from 
the (Mjnin, Mrep, Mjnax) values of the IRAC-detected sample are shown as the blue, yellow 
and red solid curves, respectively. If we assume that the stellar mass distribution of the 
blended sample is the same as the non-blended sample, we can obtain the total contribution 
from the combined non-blended and blended samples by scaling the results that we obtained 
for the non-blended sample. This scaling factor is 1.74. The dot-dashed curves in Fig. 11 
show the observed number density after this correction being applied, and are what will be 
used to compare to the models. The flattening of the observed curves at lower masses is a 
consequence of increasing incompleteness at fainter magnitudes. 

At the mass threshold of M > 1.6 x 10^°Mq, the cumulative number density derived from 
our observations is (2.6, 4.7, 10.6) xlO~^Mpc~^ based on (Mmin, Mj-ep, Mmax), respectively. 
The SPH model prediction is ~ 1.3x, 0.7x and 0.3x of the observations. The TVD model 
prediction, on the other hand, is ~ 3.7x, 2.1x and 0.9 x of the observations. This comparison 
shows that the number density of the observed high mass galaxies agrees resonably well 
with the models to within the uncertainty of the stellar mass estimates. The prediction 
from the TVD model is considerably higher than that from the SPH model, which reflects 
the uncertainty in simulations caused by the differences in the implementation methods, 
such as the mass and spatial resolutions, the simulation box sizes, and the initial conditions 
(Nagamine, priv. communications). In particular, the smaller box size of the TVD simulation 
has a significant impact at the high-mass end, as the statistics is dominated by only a small 
number of objects (only 2-3 galaxies above our choosen mass threshold). Nevertheless, the 
comparison between our observations and these simulations does show that the models of 
the hierachical paradigm are capable of producing sufficient number of high mass galaxies 
at the early epoch stage of the universe. 



5.2. Global Stellar Mass Density 

Based on the stellar masses obtained for our sample, here we estimate the global stellar 
mass density at z ^ 6. To address this question, ideally one would need a mass-limited, 
complete sample. However, the i/rs-dropout sample that we have been using is not mass- 
limited. Since it was selected from the rest-frame UV luminosities of the objects, which do 
not correlate directly with stellar masses. While the incompleteness of this sample can be 
quantified as a function of rest-frame UV luminosity, this cannot easily be transferred to the 
mass-domain. In fact, the discussion in §5.1 suggests that our sample could possibly suffer 
from significant incompleteness even at the high-mass end. Nevertheless, we can still use our 
current sample to constrain the lower limit of the global stellar mass density at z ~ 6. 
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Adding the (Afmin, M-ep, ^max) estimates of the individual objects in the IRAC-detected 
samples, the total stellar mass density locked in these galaxies is (5.8, 9.1, 31.4)xlO^M0Mpc~^. 
The contribution inferred from the IRAC-invisible samples is much smaller. The total stel- 
lar mass contributed by this sample can be approximated by multiplying the total number 
of objects to the stellar mass estimates of the median stack described in §4.2, which cor- 
responds to stellar mass density of (0.14, 0.20, 5.9)xlO^M0Mpc~^. Therefore, the sum 
of both samples is (5.9, 9.3, 37.3)xlO^M0Mpc^'^. We also correct for the incompleteness 
caused by crowding by considering the IRAC-blended sample in the same way as in §5.1 
(i.e., multiplying by a factor of 1.74). The final stellar mass density after this correction is 
(10.4,16.3,64.8) X lO^M^Mpc-^. 

This result is shown in Fig. 12, together with the estimates at lower redshifts (Dickinson 
et al. 2003). The data are shown in terms of ratios to the critical mass density, which is 
1.4xlO^^M0Mpc~^ in our adopted cosmology. We note that, however, these results have not 
yet taken into account the fact that our i775-dropout sample could still be contaminated by 
low-redshift galaxies that are not as extreme as the lEROs (see §3). Considering this effect, 
the results presented above could be reduced by ~ 23%. 

Finally, we emphasize that, while the objects in our IRAC-invisible sample make a small 
contribution to the total stellar mass density, this does not necessarily mean that the less 
massive galaxy population only makes up an insignificant fraction of the total stellar mass 
density at z ~ 6. If the MF of galaxies at z ~ 6 has a very steep slope at the low-mass 
end (see Fig. 11), the contribution from the less massive galaxies that are missing from our 
IRAC data could be much higher than that from the galaxies at the high-mass end. 



5.3. Progenitors of High-mass Galaxies: Detectability at 2; > 7 

In §4.4, we argue that the objects in the IRAC-detected sample should have much higher 
SFRs in the past. This implies that the progenitors of these objects were once much more 
UV- luminous than they are at z ^ 6. Since there have already been some efforts to push 
the current redshift boundary to 2; > 7 (e.g., Bouwens et al. 2004), it is relevant to discuss 
the prospect of detecting such early galaxies from a different perspective. 

Consider a typical object in the high-mass subset of the IRAC-detected sample that has 
Mrep = 2.0 X IO^^Mq. The stellar mass and the instantaneous SFR that its progenitor has 
at a given age depend on its SFH, as does its luminosity. Using the BC03 models that we 
adopt, we predict the apparent magnitude of its progenitor at any given redshift for a range 
of formation and observed redshifts. Fig. 13 shows two examples by assuming Zf = 7.8 (top 
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panel) and Zf = 9.0 (bottom panel), and gives its apparent magnitude in J-band if it is 
observed 50 Myr (blue) and 100 Myr (red) after the birth. We choose these two formation 
redshifts for demonstration because the age of the object when evolving to z ~ 6 would be 
~ 200-400 Myr, which is within the range of T^^p of the objects in the high-mass subset. 
In most models that we consider, the progenitor of such a very massive object should be 
detected at J ~ 24 mag. Based on the number of similar objects in the IRAC-detected 
sample, the surface density of such progenitors would be ~ 0.05 arcmin"^. While it is too 
low for the existing deep, pencil-beam surveys using the Near Infrared Camera and Multi 
Object Spectrometer onboard the HST (e.g., Thompson et al. 1999, 2005; Dickinson et 
al. 2000; Bouwens et al. 2004, 2005), such a surface density is high enough to produce a 
significant number of detections through a deep survey over a few hundred arcmin^ that 
is now feasible using the new generation of wide-field near-IR instruments at ground-based 
observatories. 

This simple argument echos the ones made about three decades ago regarding the pro- 
genitors of the local giant ellipticals, when it was suggested that such progenitors ("primeval 
galaxies") could have formed most of their stars within a short period of time and thus 
their rest-frame UV luminosities could be so high that they could be easily detected at very 
high redshifts (e.g.. Partridge & Peebles 1967; Sunyaev, Tinsley & Meier 1978; Shull & Silk 
1979). Even at that time, however, it was realized that dust extinction could easily quench 
the UV fluxes and thus those progenitors might not be as bright as expected (e.g., Kaufman 
1976). The same is also true in our discussion here. For example, a modest amount of dust 
equivalent to a reddening of E{B — V)^ 0.2 mag could easily suppress our predicted J-band 
brightness by more than 2 magnitudes, and therefore would make the detection of such pro- 
genitors aX z > 7 much more difficult. A null detection from a deep, wide-field survey would 
still be valuable, however, as it could then be used to constrain how dusty the very early 
galaxies could be and shed new light to their forming mechanism. 

There is another factor that could affect the detectability of the progenitors of high- 
mass galaxies. If the high-mass objects that we see aX z ^ 6 are the merger products of 
many much less massive progenitors, the luminosity of an individual progenitor could be 
much fainter and escape detection. In this case, a null detection would imply a high merger 
rate within the first several hundred million years of the universe, which would provide some 
constraints on the hierarchical formation theory. 
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5.4. Implication for Reionization 

As it is possible that the high-mass galaxies were born through very intense star- 
formation processes, it is interesting to study the role of their progenitors in reionizing the 
neutral hydrogen in the early universe. Specifically, here we investigate if their progenitors 
could have provided sufficient ionizing photons to sustain the reionization. To maximize their 
contribution to the ionizing photon budget, we consider an extreme - probably unphysical 
- case where all such high- mass galaxies were formed at the same time at Zy-don, which is the 
redshift when the reionization began. 

The precise value of -Zreion is still not known. Assuming instantaneous reionization with 
ionization fraction oi x^ = 1 at z < Zrdon, Kogut et al. (2003) found -Zrcion = 17± 3 based on 
the first-year WMAP results. Allowing for uncertainties in the reionization history, Kogut et 
al. obtained 11 < Zreion < 30 (at 95% confidence level). The most recent, three-year WMAP 
results have refined this estimate, but the derived value of -Zrdon is still model-dependent 
(Spergel et al. 2006). Assuming that the reionization happened instantaneously ai z > 7, 
they found the maximum-likelihood value of -Zreion = 11.3. For the purpose of demonstration, 
we use -Zreion = 9.0, 11.3 and 15.0 as examples in the discussion below. 

According to Madau, Haardt & Rees (1999; their Eqn. 27), the critical SFR den- 
sity required to reionize the universe can be calculated as p'spjiiz) ^ 0.013/g~J x [(1 -|- 
z)/6Y ^0 yr~^ Mpc~^ (assuming a Salpeter IMF). This quantity, calculated by assum- 
ing fesc = 0.1, is shown as the black, solid curves in Fig. 14 for the three assumed -Zrdon 
values. The contribution from the progenitors of the high-mass galaxies to the total SFR 
density depends on their SFH, and is less if the SFH is more prolonged. Fig. 14 shows 
the results for three toy models of exponentially declining SFHs with r of 0.1 (blue), 1.0 
(green) and 10 Myr (red). The solid, dashed and dotted curves represent the results based 
on the minimum, representative and maximum global stellar mass densities shown in Fig. 
12, respectively. These calculations demonstrate that, unless the bulk of their stars are 
simultaneously formed through an intense burst of very short duration (r < 10 Myr), the 
contribution from the progenitors of the galaxies similar to those in our sample is not suffi- 
cient to reionize the universe. Even if they indeed came into being through an unphysically 
short, simultaneous burst with r < 10 Myr, they cannot keep the universe ionized for long 
(only < 3 Myr). 

This is not to say, however, that star-forming galaxies cannot reionize the universe, 
because here we only count the progenitors of galaxies that have comparable stellar masses 
to those in our sample. If the MF has a very steep slope at the low-mass end, dwarf galaxies 
(which are below the detectability of our current observations) can contribute sufficient 
amount of ionizing photons, similar to the arguments made by Yan & Windhorst (2004a, b) 



-20- 

for the case at z ~ 6. 

6. Summary 

In this paper, we use the full-epoch IRAC observations of the GOODS Spitzer Legacy 
Program to study the stellar masses of galaxies aX z ^ 6. In total, 53 iyys-dropouts (seven of 
them with spectroscopic confirmation) selected in the entire GOODS fields (~ 330 arcmin^) 
are securely detected by the IRAC. We derive the stellar masses for all these objects, using 
three different mass estimators: Mmjn, Mj-ep, and Mmax- We argue that the true mass of 
an object should be between M^i^i and Mmax, and use Mrep as its "representative mass". 
We further refer to the age to which Mrep corresponds as its "representative age", T^ep- 
We find that the Mrep values range from 0.09-7.0 xlO^^M©, and the Trep values range from 
50-400 Myr, both are consistent with the results of our earlier work based on a smaller 
sample in the HUDF (Paper I). We also study 79 iyys-dropouts that are invisible from the 
IRAC observations, and find that they are generally one order of magnitude less massive 
than those detected by IRAC. Their ACS-to-IRAC colors are much bluer than those of the 
IRAC-detected ones, indicating that their luminosities are dominated by young population 
of unreddened stars with ages < 40 Myr. We discuss various sources of uncertainty in the 
mass estimates, and find that most of our mass estimates are relatively robust. 

We discuss a number of implications of our results. As in Paper I, we find that the 
existence of the most massive galaxies in our much enlarged sample can still be explained 
by at least one set of N-body simulations of the hierarchical paradigm. We derive a lower 
limit to the global stellar mass density at z ^ 6, and conclude that at least 0.2-1.1% of 
the total stellar mass in the local universe has been locked in stars hj z ^ 6. As the most 
massive galaxies in our sample should have much larger SFRs in the past, we discuss the 
prospect of detecting their progenitors at 2; > 7. A near-IR survey covering a few hundred 
arcmin^ to a depth of AB ~ 24.0 mag can provide valuable constraints on the properties of 
such progenitors even in the case of a null detection. We also discuss how the progenitors 
of galaxies comparable to those in our sample could contribute in reionizing the hydrogen in 
the early universe, and conclude that such progenitors alone are not sufficient to sustain the 
reionization of the universe. 

We thank the other members of the GOODS team who have contributed to the success 
of the observations and data analysis. We also thank our referee, Rodger Thompson, for 
his prompt and insightful referee report. We are grateful to Kentaro Nagamine for his up- 
to-date MP. We wish to thank Ivo Labbe for useful discussion. Support for this work, part 



-21- 



of the Spitzer Space Telescope Legacy Science Program, was provided by NASA through 
Contract Number 1224666 issued by the Jet Propulsion Laboratory, Cahfornia Institute of 
Technology under NASA contract 1407. HY acknowledges the support from the NASA grant 
HST-GO-09780.03. 



REFERENCES 

Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393 

Brinchmann, J. & Elhs, R. S. 2000, ApJ, 536, L77 

Bouwens, R. J., et al. 2004, ApJ, 616, L79 

Bouwens, R. J., Illingworth, G. D., Thompson, R. L, Franx, M. 2005, ApJ, 624, L5 

Bruzual, A. G. & Chariot, S. 2003, MNRAS, 344, 1000 (BC03) 

Calzetti, D., Armus, L., Bohlin, R. C, Kinney, A. L., Koomneef, J., & Storchi-Bergmann, 
T. 2000, ApJ, 533, 682 

Chabrier, G. 2003, PASP, 115, 763 

Chary, R., Stern, D. & Eisenhardt, P. R. M. 2005, ApJ, 635, L5 

Chen, H.-W. & Marzke, R. O. 2004, ApJ, 615, 603 

Cohen, J. G. 2001, AJ, 121, 2895 

Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168 

Dickinsion, M., et al. 2000, ApJ, 531, 624 

Dickinson, M., Paopvich, C, Ferguson, H. C, & Budavari, T. 2003, ApJ, 587, 25 

Dickinson, M., et al. 2004, ApJ, 600, L99 

Egami, E., et al. 2005, ApJ, 618, L5 

Eyles, L. P., Bunker, A. J., Stanway, E. R., Lacy, M., Ellis, R. S., & Doherty, M. 2005, 
MNRAS, 364, 443 

Fazio, G. G., et al. 2004, ApJS, 154, 10 

Giavalisco, M., et al. 2004a, ApJ, 600, L93 



-22- 

Giavalisco, M., et al. 2004b, ApJ, 600, L103 

Kaufman, M. 1976, Ap&SS, 40, 369 

Kroupa, P. 2001, MNRAS, 322, 231 

Kogut, A., et al. 2003, ApJS, 148, 161 

Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648 

Madau, P., Pozzetti, L. & Dickinson, M. 1998, ApJ, 498, 106 

Mobasher, B., et al. 2005, ApJ, 635, 832 

Nagamine, K., Cen, R., Hernquist, L., Ostriker, J. P., & Springel, V. 2004, ApJ, 610, 45 

Night, C, Nagamine, K., Springel, V., & Hernquist, L. 2006, MNRAS, 366, 705 

Papovich, C, Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620 

Partridge, R. B & Peebles, P. J. E. 1967, ApJ, 147, 868 

Rudnick, G., et al. 2003, ApJ, 599, 847 

Salpeter, E. E. 1955, ApJ, 121, 161 

Shull, J. M. & Silk, J. 1979, ApJ, 234, 427 

Sunyaev, R. A., Tinsley, B. M., & Meier, D. L. 1978, ComAp, 7, 183 

Spergel, D. N., et al. 2003, ApJS, 148, 175 

Spergel, D. N., et al. 2006, submitted to ApJ (astro-ph/0603449) 

Stanway, E. R., et al. 2004a, ApJ, 604, L13 

Stanway, E. R., Bunker, A. J., McMalion, R. G., Ellis, R. S., Treu, T., & McCarthy, P. J. 
2004b, ApJ, 607, 704 

Thompson, R. I., Storrie-Lombardi, L. J., Weymann, R. J., Rieke, M. J., Schneider, G., 
Stobie, E., Lytle, D. 1999, AJ, 117, 17 

Thompson, R. I., et al. 2005, AJ, 130, 1 

Vanzella, E., et al. 2006, submitted to A&A (astro-ph/0601367) 

Werner, M. W., et al. 2004, ApJS, 154, 1 



-23- 

Yan, H., et al. 2004, ApJ, 616, 63 
Yan, H., et al. 2005, ApJ, 634, 109 (Paper I) 
Yan, H. & Windhorst, R. A. 2004a, ApJ, 600, LI 
Yan, H. & Windhorst, R. A. 2004b, ApJ, 612, L93 



This preprint was prepared with the AAS I^TJtjX macros v5.0. 



-24- 



40 - 



30 - 



I 20 h 



10 - 







1 — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 



□ All 

IRAC-detected 
IRAC-invisible 



24 



f^ R^ 



25 




v^ ri 



29 



'850 



Fig. 1. — The distribution of ^sso magnitudes of the irrs-dropouts in the GOODS fields. 
The shaded blue and red histograms are for the IRAC-detected and IRAC-invisible samples, 
respectively, while the black, unfilled histograms are for all objects (including the blended 
sample). 
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Fig. 2. — This figure demonstrates the evolution of niass-to-light-ratio {M/L; shown here 
in terms of M/Ly) for the models considered in our study (see §4). The red curve at the 
top corresponds to a SSP, the green curves correspond to models with r ranging from 10 to 
90 Myr, while the blue curves correspond to models with r from 0.1 to 1.0 Gyr. At the age 
of ~ 1 Gyr, which is the maximum age allowed by our adopted cosmology, a SSP reaches 
the largest possible M/L among all the models, and thus gives the highest possible stellar 
mass for a given luminosity. 
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Fig. 3. — The evolution of {z^^q — m^Q^m) color for the models considered in our study (see 
§4). The curves are color-coded in the same way as in Fig. 2. The {zs5o — fn^Sfim) color can 
be used as an age indicator. 
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Fig. 4. — The top panel displays the histograms of the stellar mass estimates of the galaxies 
in the IRAC-detected samples, using three different estimators as described in §4.1. The 
vertical dashed line in each figure represents the corresponding mass estimate of a typical 
object in the IRAC-invisible sample (see §4.2). The bottom panel shows the histogram of 
the upper limits of the stellar masses of the galaxies in the IRAC-invisible sample, calculated 
using their flux density upper limits (2 a) in the 3.6/im-channel (see §4.2). 
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Fig. 5. — The (2:850 — mssfim) colors of the IRAC-detected objects are shown as red sohd 
squares, while the upper limits of these colors of the IRAC-invisible objects are shown as 
downward arrows. The blue hexagon represents the median color of the IRAC-invisible 
objects, which is calculated by using the median zs5o magnitude of the sample and the 
m^Q^rn magnitude of the median stack. The IRAC-invisible objects are significantly bluer 
than the IRAC-detected ones, indicating that the young populations are playing a more 
important role in them than in the IRAC-detected objects. 
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Fig. 6. — The 3.6/im median stack of the objects in the IRAC-invisible sample is displayed 
in the left panel, where there is a visible source in the middle (indicated by the circle). For 
comparison, the median stacks of the same number of random positions is shown in the right 
panel and shows no detected source in the middle circle. 
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Fig. 7. — The systematic uncertainties of the stellar mass estimates caused by fixing the 
redshifts of our models at 2; = 6. The evaluations of the three mass estimates are repeated 
for every IRAC-detected object by using models a.t z = 5.5 (asterisks) and z = 6.5 (squares). 
The systematics are expressed in terms of relative differences of Mmin, Mj-^p and Mmax in the 
top, middle and bottom panels, respectively. We define AM^ = M — M , where M is the 
value obtained with the original models at z = 6, and M is the value obtained with the 
new models. Thus a positive data point means that our original evaluation overestimates 
the quantity if the object is actually at a redshift different from z = 6, and a negative data 
point means the opposite. If our objects are actually all at 2; = 5.5, our original results 
underestimate Mmin and M^^p only by a few percent on average, but overestimate M^ax by 
21%. If our objects are actually all at z = 6.5, our original results overestimate both Mmin 
and Mj-cp by ~ 25%, and underestimate Mmax by 27%. 
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Fig. 8. — Similar to Fig. 7, but for the systematic uncertainties caused by using models of 
solar metallicity only. The mass estimates are recalculated by using models with metallicity 
of 1/200 of solar, which is the lowest abundance available for the BC03 models. The relative 
difference of mass is defined in a similar way to the quantity shown in Fig. 8. If all our 
objects actually have such a low metallicity, our original results underestimate Mmin and 
M^ep by ~ 70% and 32% on average, respectively, and overestimate Mmax by 32%. 
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Fig. 9. — Similar to Fig. 7 and 8, but for the systematic uncertainties caused by using zero- 
reddening models. The mass estimates are recalculated by using models with E{B — V) = 0.2 
mag, following the reddening law of Calzetti et al. (2000). The relative difference of mass is 
defined in a similar way as before. If all our objects actually suffer from dust reddening and 
obscuration of this amount, our original results overestimate Mmin and Mj-ep by ~ 13% and 
30%, respectively, and underestimate Mmax by ~ a factor of 3.35. 
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Fig. 10. — The specific SFRs of tlie IRAC-detected objects (filled triangles) as a function 
of Mrep. For completeness, the specific SFR of the IRAC-invisible objects, calculated using 
their median zgso magnitudes and the median 3.6/im stack, is also shown (the asterisk). 
The two solid lines shows the specific SFRs of objects with constant SFRs of 3.0 (left) and 
30 Mq yr~^ (right). The label to the right shows the time periods through which the galaxies 
can assemble their total masses if their past SFRs are the same as the current values. The 
horizontal dashed line indicates the age of the universe at z ^ 6. Objects below this line 
must have acquired their masses at a much higher SFR than the current values. 
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Fig. 11. — The comparison between the observed number density inferred from our sample 
and the mass functions predicted by a set of N-body hydrodynamic simulations in a ACDM 
universe from Nagamine et al. (2004) and Night et al. (2006). Both are presented as 
cumulative number densities. The simulations are shown as the black solid (SPH) and 
dashed (TVD) curves, while the observed number densities based on (Mmin, Mj-ep, Mmax) 
are shown as the blue, yellow and red curves, respectively. The solid, color-coded curves 
are the observed values based only on the non-blended, IRAC-detected sample, while the 
dot-dashed, color-coded curves are those after taking the blended sample into account. The 
asterisk to the left is the observed cumulative density when assuming all the IRAC-invisible 
objects (after correction for blending) have the M^p value as derived using the median stack 
(see §4.2). The dotted vertical line at (M = 1.6 x lO^^M©) indicates the mass threshold 
above which our IRAC data (but not necessarily the irys-dropout sample itself) are complete. 
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Fig. 12. — The lower limit to the global stellar mass density at 2; ~ 6 based on our study is 
shown as the rectangle in this figure. The open star inside this rectangle shows the result 
based on the Mrep values, while the lower and upper boundaries are set based on the Mmin 
and Mmax values, respectively. The upward arrow indicates that these results are lower 
limits. The results at lower redshifts are from Cole et al. (2000), Brinchmann & Ellis (2000), 
Cohen (2001), Dickinson et al. (2003) and Rudnick et al. (2003). The Y-axis labels to right 
are in unit of Mq/Mpc^, while the ones to the left are in terms of ratio to the critical mass 
density, which is 1.4 x IO^^MqMpc"'^ for our adopted cosmology. 
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Fig. 13. — The predicted J-band brightness of the progenitor of a typical high-mass object 
(M = 2 X IO^^Mq) in the IRAC-detected sample. The brightness depends on how quickly 
the progenitor has been forming stars (characterized by the timescale r of the exponential 
SFH) and when it is observed after its birth. Two cases, 50 and 100 Myr after the birth, are 
shown as the blue and red curves, respectively. With a formation redshift of zj = 7.8 (top 
panel), these two ages correspond to being observed at z ~ 7.4 and 7.0, respectively, while 
with a formation redshift of zj = 9.0 (bottom panel), they correspond to being observed at 
z ^ 8.4 and 8.0, respectively. These two formation redshifts are chosen for demonstration 
because a galaxy evolving to 2; ~ 6 would have an age of 200-400 Myr, which spans the T^ep 
range of the M > IO^^Mq objects in the IRAC-detected sample. The calculation is done for 
r = 10 Myr to 1 Gyr based on the BC03 models, assuming no dust extinction. An modest 
amount of dust corresponding to E{B — V) ^ 0.2 mag would lower the brightness by ~ 2 
mag. 
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Fig. 14. — The contribution to the reionizing photon background by primeval galaxies, 
inferred from the progenitors of the objects in our sample. We assume an extreme case where 
the universe started to give birth to all the progenitors of the galaxies in our sample at the 
same time at ^^reion, and show the contribution of these progenitors to the reionizing photon 
background, expressed in terms of SFR density, as a function of time after the birth. This 
contribution is very sensitive to the SFH of the galaxies. Here we show the results for three 
exponentially declining SFHs with r of 0.1, 1.0 and 10 Myr as blue, green and red curves, 
respectively. The solid, dashed and dotted curves represent the results that correspond 
to the minimum, representative and maximum global stellar mass densities shown in Fig. 
12, respectively. The black curves, corresponding to Zreion = 15.0, 11.3 and 9.0, are the 
critical SFR densities needed to fully ionize the neutral hydrogen in the universe, calculated 
according to Madau, Haardt & Rees (1999). See §5.2 for detailed interpretation. 



